/* 
   indexg_gc - Build an gc content index for the whole genome

   Copyright (C) 2007 

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software Foundation,
   Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

   #include <termios.h>
   #include <grp.h>
   #include <pwd.h>
*/

#include <stdio.h>
#include <sys/types.h>
#include <getopt.h>
#include "system.h"
#include "common.h"

/* include my fastaio lib */
#include "fastaIO.h"

#define EXIT_FAILURE 1

static void usage (int status);

/* The name the program was run with, stripped of any leading path. */
char *program_name;

/* getopt_long return codes */
enum {DUMMY_CODE=129
};

/* Option flags and variables */

static char * ifile = NULL;		/* --input */
static char *oname = "stdout";		/* --output */
static UBYTE2 window = 200;		/* --window */
static UBYTE2 step = 10;		/* --step */
static FILE * fhifile; 			/* input file handle */
static FILE *ofile;

static struct option const long_options[] =
  {
	{"input", required_argument, 0, 'i'},
	{"window", required_argument, 0, 'w'},
	{"step", required_argument, 0, 's'},
	{"help", no_argument, 0, 'h'},
	{"version", no_argument, 0, 'V'},
	{NULL, 0, NULL, 0}
  };


static int decode_switches (int argc, char **argv);

/* buffer read  */
char *
buffer_push (FILE * fp, char * interval)
{
  long i;
  long l;
  l = strlen(interval);
  for (i=0;i<l-step;i++) {
	interval[i] = interval[i+step];
  }
  for (i=l-step;i<l;i++) {
	interval[i]=fgetc(fp);
  }
  return interval;
}

BYTE
cal_gc ( char * seq)
{
  BYTE gc;
  long gcsum;
  int i;

  gcsum = 0;
  for (i=0;i<window;i++) {
	if (seq[i] == 'N')
	  return -1;
	if (seq[i] == 'C'||seq[i] == 'G'||seq[i] == 'c'||seq[i]=='g')
	  gcsum ++;
  }
  gc = gcsum*100/window;
  return (BYTE)(gc);
}



int
main (int argc, char **argv)
{
  int i;
  long char_num;
  UBYTE4 len;
  char *seq, *text, *id;
  BYTE gc;
  FastaRecord * records;
  FastaRecord * cur_record;
  FastaRecord * tsub_record;


  program_name = argv[0];

  i = decode_switches (argc, argv);

  if ( ! fhifile )
    usage(0);

  /* do the work */

  /* read all sequences */
  fseek(fhifile,0,SEEK_END);
  char_num = ftell(fhifile);
  rewind(fhifile);
  text = (char*) malloc(sizeof(char)*(char_num+1));
  fread(text,sizeof(char),char_num,fhifile);
  fclose(fhifile);
  /* parse the whole text, produce a chain of fasta records */
  records = parse_text (text);
  free(text);

  cur_record = records;

  while (cur_record)
	{
	  up_record(cur_record);
	  len = cur_record->fasta->length;
	  id = cur_record->fasta->id;

	  fwrite (id,sizeof(char),10,ofile);
	  fwrite (&window,sizeof(UBYTE2),1,ofile);
	  fwrite (&step,sizeof(UBYTE2),1,ofile);
	  fwrite (&len,sizeof(UBYTE4),1,ofile);
	  seq = cur_record->fasta->seq;
	  for (i=0;i<((len-window)/step+1);i++) {
		gc = cal_gc(seq);
		fwrite (&(gc),sizeof(BYTE),1,ofile);
		seq+=step;
	  }
	  cur_record = cur_record->next;
	}
  free_records(records);
  fclose(ofile);

  exit (0);
}

/* Set all the option flags according to the switches specified.
   Return the index of the first non-option argument.  */

static int
decode_switches (int argc, char **argv)
{
  int c;

  ofile = stdout;

  while ((c = getopt_long (argc, argv, 
						   "i:"	/* input */
						   "o:" /* output */
						   "w:"	/* window */
						   "s:"	/* step */
						   "h"	/* help */
						   "V",	/* version */
						   long_options, (int *) 0)) != EOF)
    {
      switch (c)
		{
		case 'i':
		  ifile = strdup(optarg);
		  fhifile = fopen(ifile,"r");
		  if (!fhifile) {
			fprintf (stderr, "cannot open `%s' for reading\n", ifile);
			exit(1);
		  }
		  break;
		case 'o':		/* --output */
		  oname = (char *)strdup(optarg);
		  ofile = fopen(oname, "wb");
		  if (!ofile)
			{
			  fprintf(stderr, "cannot open `%s' for writing\n", oname);
			  exit(1);
			}
		  break;
		case 'w':
		  window = strtod(optarg,NULL);
		  break;
		case 's':
		  step = strtod(optarg,NULL);
		  break;
		case 'V':
		  printf ("indexg_gc %s\n", VERSION);
		  exit (0);
		  
		case 'h':
		  usage (0);

		default:
		  usage (EXIT_FAILURE);
		}
    }
  
  return optind;
}


static void
usage (int status)
{
  printf (_("%s - \
Calculate GC content for every intervals along the whole genome\n"), program_name);
  printf (_("Usage: %s <-i fasta file> [-o output file] [OPTION]... [FILE]...\n"), program_name);
  printf (_("\
Options:\n\
  -i, --input                fasta file to read\n\
  -o, --output               output binary file, default: stdout\n\
  -w, --window               window size, default: 200\n\
  -s, --step                 step size, default: 10\n\
  -h, --help                 display this help and exit\n\
  -V, --version              output version information and exit\n\
\n\
Note:\n\
The Output Binary File is divided in blocks, the format of every block is:\n\
Bytes\tContent\n\
1-10\tChromosome Name\n\
11-12\tWindow size\n\
13-14\tStep size\n\
15-18\tLength of Chromosome\n\
19-end of block\tGC for every Position till the end(-1 is for NA)\n\
\n\
Author:\n\
  Tao Liu: taoliu@jimmy.harvard.edu\n\
"));
  exit (status);
}
